#### Produce figures describing municipality-level coefficients

### Prepare environment

setwd("~/research/coffee-dataverse/src")

do.origspec <- T # For determining priors
source("load.R", encoding="iso-8859-1")

library(rstan)
library(dplyr)

output.suffix <- "-eq24-orig"

## Classify according to portion arabica
df.mean <- df %>% group_by(Munic�pio) %>% summarize(arabica.portion=mean(ifelse(is.na(arabica.produce), ifelse(is.na(robusta.produce), NA, 0), arabica.produce) / total.produce, na.rm=T), total.average=mean(total.produce, na.rm=T), tmax=mean((tmax01 + tmax02 + tmax03 + tmax04 + tmax05 + tmax06 + tmax07 + tmax08 + tmax09 + tmax10 + tmax11 + tmax12) / 12, na.rm=T), tmin=mean((tmin01 + tmin02 + tmin03 + tmin04 + tmin05 + tmin06 + tmin07 + tmin08 + tmin09 + tmin10 + tmin11 + tmin12) / 12, na.rm=T))

## Collect all municipality-level results

param <- get.fit.params(output.suffix)

betas <- list()
output.suffixes <- c("-eq24-orig", "-arabica-eq24-orig", "-robusta-eq24-orig")
for (output.suffix in output.suffixes) {
    load(paste0("../data/combined-betagamma", output.suffix, ".RData"))
    means <- apply(la$beta, 2:3, mean)
    colnames(means) <- param

    muniinfo <- read.csv(paste0("../data/combined-muniinfo", output.suffix, ".csv"))
    allinfo <- cbind(muniinfo, means)

    allinfo2 <- allinfo %>% left_join(df.mean, by=c('region'='Munic�pio'))

    betas[[output.suffix]] <- allinfo2
}

## Display graphs for each parameter

library(ggplot2)

labels <- list('death_coeff1'="Die-off EDD effect (per 1000 EDDs)", 'autoreg_planting'="Planting autoreg.", 'autoreg_price'="Nerlove autoreg.", 'cost_harvest'="Cost of harvest", 'cost_to_planting'="Planting price resp.", 'scale_trend'="Exogenous trend", 'yield_sigma'="Yield std. dev.", 'yield_uncertainty'="Range of yields", 'bearing_sigma'="Bearing std. dev.", 'harvest_sigma'="Harvest std. dev.", 'production_sigma'="Production std. dev.")
for (kk in 1:length(weatherlabels)) {
    labels[[weatherpreds[kk]]] <- weatherlabels[kk]
    labels[[paste0('yield_coeff', kk)]] <- weatherlabels[kk]
}
labels[[paste0('yield_coeff', which(weatherpreds == 'gdd1000'))]] <- "Yield GDD effect (per 1000 GDDs)"
labels[[paste0('yield_coeff', which(weatherpreds == 'edd1000'))]] <- "Yield EDD effect (per 1000 EDDs)"

for (pp in 1:length(param)) {
    pdf.all <- betas[['-eq24-orig']][, c('region', 'tmin', 'tmax', 'arabica.portion', param[pp])]
    pdf.arabica <- betas[['-arabica-eq24-orig']][, c('region', param[pp])]
    pdf.robusta <- betas[['-robusta-eq24-orig']][, c('region', param[pp])]

    pdf <- pdf.all %>% left_join(pdf.arabica, by='region', suffix=c('', '.arabica')) %>%
        left_join(pdf.robusta, by='region', suffix=c('', '.robusta'))

    names(pdf)[5:7] <- c('param', 'param.arabica', 'param.robusta')

    ggplot(pdf, aes(tmax - 273.15, param, colour=arabica.portion)) +
        geom_point() + geom_smooth(colour='black') + theme_bw() +
        scale_x_continuous(expand=c(0, 0)) + xlab("Average Maximum Temperature") +
        ylab(labels[[param[pp]]]) + theme(legend.position="bottom") +
        scale_colour_gradient2(name="Portion Arabica:", low='#d7191c', high='#1a9641', mid='#bfbf7f',
                               midpoint=.5, labels=scales::percent)
    ggsave(paste0("../paper/bymuni/", param[pp], '-vs-tmax.pdf'), width=6.5, height=6.5)

    ggplot(pdf, aes(tmin - 273.15, param, colour=arabica.portion)) +
        geom_point() + geom_smooth(colour='black') + theme_bw() +
        scale_x_continuous(expand=c(0, 0)) + xlab("Average Minimum Temperature") +
        ylab(labels[[param[pp]]]) + theme(legend.position="bottom") +
        scale_colour_gradient2(name="Portion Arabica:", low='#d7191c', high='#1a9641', mid='#bfbf7f',
                               midpoint=.5, labels=scales::percent)
    ggsave(paste0("../paper/bymuni/", param[pp], '-vs-tmin.pdf'), width=6.5, height=6.5)
}
